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Motivated by recent experiments by Yuse and Sano (Nature, 362, 329 
(1993)), we propose a discrete model of linear springs for studying fracture in 
thin and elastically isotropic brittle films. The method enables us to draw a 
map of the stresses in the material. Cracks generated by the model, imposing 
fSj ' a moving thermal gradient in the material, can branch or wiggle depending 

on the driving parameters. The results may be used to compare with other 

recent theoretical work, or to design future experiments. 
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I. INTRODUCTION. 



The study of brittle fracture is of particular technological interest for understanding 
the properties of ceramic materials. By ceramic we understand covalent ionic materials, 
including also glasses, policrystaline aggregates, minerals and even composites |]]. Studies 
in the statistical physics and engineering aspects of the problem 0|| have been carried out, 
and though some proposals have been made, a comprehensive and general understanding 
of the mechanisms leading to fracture is still missing. The fundamental reason for this 
apparent lack of progress may be that experimental work in the study of these processes is 
very difficult because of the great quantity of parameters that seem to have impact in the 
resulting shape or velocity of the crack. Recently Yuse and Sano carried out an experiment 
H in which they could control the fracture pattern as function of the external imposed 
stresses. It is sketched if Fig. 1. A thin rectangular plate of glass with an initial notch is 
pulled at a velocity v inside a region of width 2£, where there is a change in the temperature 
from Th to T c . As function of increasing v, and also AT=T^ — T c , several phases may be 
found: first, for low values of the control parameters, no fracture is propagated. If the 
control parameter, i.e. v, is increased, a straight crack is generated. If it is again increased, 
a wiggling pattern similar to that of Fig. 1 appears. And, finally, for high enough values of 
v branching occurs, and more than one crack propagates in the system. 

This experiment has inspired directly analytical |5||| as well as numerical |7],[8| work. 
Recent numerical simulations |S| also show this oscillatory instability when the system is 
driven, not by a thermal gradient, but by an increasing strain in the borders. 

The purpose of this paper is to introduce a discrete model that has the right long wave 
behavior of the continuum elasticity theory, and to check its predictions of wiggling cracks 
in the case of a thermally driven fracture. Here I will thus try to present the tool. Quan- 
titative comparisons against analytical models and experiment will be the subject of future 
publications. 

The basic physical ideas that inspire the model developed here are that brittle materials 
are elastic, in the strict sense that displacements are proportional to applied stresses, until 
they fail. Also, since the velocity of the crack propagation is imposed by the motion of 
the temperature change, and it is much smaller than the sound velocity in the material, 
cuasistatic calculations are good candidates to give a right description of the problem. 

In section II I will describe the model, expose the fracture processes, and see how the 
model reproduces qualitatively the thermally driven fracture scenario. In section III I will 
show how the stress profiles in the system can be calculated. Section IV will be dedicated to 
study the dependence of the different parameters when we want to produce larger systems 
(or better discretizations). Finally, in section V, we will present some conclusions and paths 
for future works and applications of the model. 



II. MODEL AND SIMULATIONS. 
A. Discretization. 

The continuum elastic description of a two dimensional material that is in the presence 



of a temperature field, T, is given by the free energy [10 
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F e = ^(A + fi){e xx + Eyy) 2 + fi{2e 2 xy + -(e xx - e yy ) 2 ) - (A + fi)aT(e xx + e yy ), (1) 

where = (diUj + djUi)/2 are the components of the strain tensor (v,i are the local defor- 
mations of the material in the % direction), A and fi the Lame coefficients, and a the thermal 
expansion. 

As simple microscopic model we can consider a system of hookean springs described by 
the free energy 

= \Yj k a t( u i - u *) ■ - E k a w i ( u i - u *) ■ % ( 2 ) 

ij ij 

where Uj is the displacement of node i from its equilibrium position, is the unitary vector 
that points from node i to node j, kij is the force constant of the spring that joins zth and 
jth nodes, and W{ is the scalar field that includes temperature, compression modulus, and 
thermal expansion at node i [11]. If the nodes are placed in a triangular lattice, and links 
are made all with the same force constant and only between nearest neighbors, the long 
wavelength behavior of Eq. |2| is the same as that described by the continuum Eq. [1], when 
the temperature is constant in the material, and with A=/z. 

Let us introduce a reference system in terms of the axis of the triangular lattice, making 
the following changes of variables: u = (2/y / 3)x and v = y + (l/y/S)x for the coordinates, 
and u u = (v / 3/2)m x — (l/2)u y , u v = u y , and u w = u u + u v for the three vector displacements 
in the triangular lattice f!2]| . We will also take the vectors fy as vectors in the triangular 
lattice, making the equations of motion linear in the displacements. 

The free energy of the model can be written as 

= 9 {(. d uU u ) 2 + {d v u v ) 2 + {d w u w ) 2 ^j - 2 W n {d u u u + d v u v + d w u w ) , (3) 

where n now sweeps all the nodes in the system, and the derivatives are evaluated as differ- 
ences between neighbor nodes when there are links between them (k^ takes values (1) for 
missing (present) springs). 



B. Temperature. 

If the temperature drop is being displaced at a velocity v in the y direction, the thermal 
diffusion equation 

accepts solutions that depend only on ip = y — v t in the form: 

rrf ,x T h -T c T h + T c cosh v^ - exp vi/j 

TW = + —2 sinh< ' (5) 

where v= v/kis the inverse of the temperature length, and 2£ is the distance between the 
hot and the cold temperature baths. Function T is also sketched in Fig. 1. Notice that its 
derivative is not continuous at T = T c , where the fastest changes in T are found, and, for 
that reason, larger stresses are located. 
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C. Initial conditions and relaxation. 



To perform a simulation the nodes are placed in a rombus of dimensions N u x N v , in 
units of the node-node distance. Then the corners are cut to make a rectangular sample of 
dimensions (v3/2)iV u x N v -N u /2. A small notch cutting bonds at the edge simulate the 
experimental initial conditions. Now the temperature drop is introduced by applying W, 
following Eq. || with y = vt far bellow the system. The set of nodes displacements, {u}, is 
then relaxed by the Conjugate Gradients Method |]I3| until forces (elastic and thermal) at 
each node are balanced. 



D. Fracture. 



Once the system is relaxed, springs whose stress (including difference in displacements 
and average W in the connected nodes) is more than some threshold value, ath, are allowed 
to break. Since pure deterministic breakdown can give rise to problems in the decision of 
which spring to break |?J, we will make use of probabilistic breakdown following Louis and 



Guinea model |14|], and considering four neighbors of the already broken springs ||15|| . In this 
model each spring has a probability of breaking proportional to a power 77 of its stress. If 
77 is large enough (I will take in this work 77=10), the method is, for all practical purposes, 
deterministic, except in those cases close to branching where there are two or several paths 
to which fracture may develop. If after the relaxation there are no springs with stresses 
larger than the threshold, the temperature field is displaced in the positive y-direction, y 
— > yo+Ay, and the system is relaxed again. If breakdown occurs, the system is also relaxed 
without the broken bond and stresses recalculated without moving the temperature profile. 
The simulation stops when the crack gets close to the upper limit of the system. 



E. Results. 

As discussed in previously developed discrete models J7|, the different shapes of cracks 
do not depend on the values of Wh and W c , but on its difference, AW. In particular, the 
model presented here, as compression or dilation generate the same stressed state, depends 
only on the absolute value of AW. The units of AW are the same as a t h, and properties 
are only dependent on AW/a t h- We will thus set a t h = 1, and discuss its possible changes 
in section IV. 

Several simulations with £ = 6, and system size 30x100 are shown in Fig. |2|. The model 
reproduces qualitatively well the phenomenology of thermally driven fracture described in 
section I. In Fig. 2 can be found several interesting configurations, ranging from branching 
in the middle of the plate, some time after the crack seems to be a wiggling one, to situations 
in which wiggling and straight cracks run parallel. In fact, it seems that, with enough CPU 
time, besides the line transitions from straight to wiggling, and from wiggling to branching, 
a line transition can be defined for any branching level, or even for transitions from two 
straight cracks to one wiggling and one straight. 

This model has not the invariance reported in Ref. @ with respect of width-v, or, in other 
words, that the effective thermal lenght, sets the units of the width of the system. This is not 
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the case here because we have considered another parameter, that is, the distance between 
the two applied temperature baths, 2£. As a consecuence, the temperature as function of 
distance is an exponential and not an hyperbolic tangent. 



III. FIELD OF STRESSES. 



Stresses are the derivatives of the free energy with respect to the strains, = dT s jdui^ 

3 

&xx = ^ (d u u u + d w u w ) - 3W, (6a) 
\/3 

Cxy = — (d w u w - d u u u ) , (6b) 



1 

a yy = ^ (duU u + d w u w ) + d v u v - 3W. (6c) 



The same equations, but for a constant factor, would have been obtained if instead of 
Eq. ^ Eq. [I] had been used. This is a result that enhances the goodities of the triangular 
symmetry in the description of homogeneous media. Eqs. 6 are evaluated at each node of 
the system, considering the partial derivatives as finite differences of the displacements of 
the neighbor connected nodes. 

To illustrate the capabilities of the method, two different stress conditions that lead to a 
straight crack are shown in Fig. [| Both system sizes are 30x100. The shading correspond 
to dividing the set of stresses in 10 greys, and putting the same amount of points in each 
color. This is done to see more clearly the symmetries in each component of the stress 
tensor, and the differences between the two drivings. The frames are taken when the crack 
tip is in the middle of the plate. The stress driven plate is simulated by applying forces at 
each boundary node such that a xx = ctq, and a xy = a yy = 0, when there are no springs cut 
from the plate. Temperature now is constant (zero). For the stress driven plate, the stress 
of the crack tip grows as crack evolves, where probably the cuasistatic description used here 
fails due to the acceleration of fracture. Nevertheless, qualitatively, the stress profiles shown 
in Fig. ^ must remain the same in a dynamic description of the problem. 

As seen from Fig. [| the symmetries in the stress tensor for the thermally driven fracture 
are more complicated than for the stress driven one. Symmetries of different stresses that 



can be also found by analytically solving the equations for the Airy function |L6| , obtaining 
results consistent with the simulations. 

To take a closer look to the crack tip divergence, sections along the middle of the plate 
are drawn in Fig. 4, where a xy is not plotted, as is, by symmetry, zero, or, numerically, very 
close to zero due to effects of the small lattices considered. Unfortunately the divergence of 
the stresses near the crack tip consist only of a couple of data points that are not enough to 
determine the exponents of the divergence. 



IV. SIZE DEPENDENCE. 

One of the objections that may be argued about last section is that the size of the cracks 
shown was too small to get any detail of questions of fundamental importance like stress 
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divergence at the tips, or even the location of the tip itself. As we are treating a model that 
is but a mere discretization of the elastic equations, if all the lengths of a given simulation 
(system size, £, and It = 1/v), where scaled by a factor, the new simulation would reproduce, 
at larger scale, the patterns obtained in the small scale one (and much less CPU consuming 
one). But be forgot to rescale the stress threshold ath, that has really units of stress divided 
by length. 

To guess this scaling we will follow an argument also used in effective medium theories 
for elastic percolation fllTjl . We can apply a force of modulus /, pulling two nodes of the 



triangular lattice apart. It would produce a displacement proportional to /, and the pro- 
portionality constant is the effective spring constant in the system, 1/a*, due to all possible 



connections of the two nodes through the infinite triangular lattice ||18|| . As a t h takes in 
account only the direct connection between neighbor nodes, the contribution of the medium 
must be subtracted. The expression for the threshold stress, when the lengths of the system 
are multiplied by n, is given, in this approximation, by 



where a* is calculated putting the force / n nodes apart in one of the directions of the 
triangular lattice. 

Eq. |7| turns out to be only a good approximation (really an upper bound) to the value 
that rescales the system to find the same crack patterns reproduced. Fig. 5 shows such 
rescaling for one of the cracks shown in Fig. 2. The shape (wavelenth and amplitude) are 
the same for the three discretizations, but the threshold stresses are not exactly those given 
by Eq. [7| (changes in ath are translated to changes in AW since in the simulations the 
convention is that a th = 1). For size 60x200 (90x300), a th = 1.452(1.8) in Fig. 5, and from 
Eq. [7] we would have that a t h — 1.475(1.952), values that lead to a straight crack in the 
simulations. 



Thus, the theoretical values for ath obtained by this method are not the right ones ||19 
but provide a good estimation that can save time when we want to have larger systems 
and we do not want to sweep the hole parameter space involving large simulations. As a 
reference, the 90x300 crack presented here took 5 days cpu time in an alpha machine. 



V. CONCLUSIONS. 

As conclusions, we have seen that hookean springs placed in a triangular lattice, and 
coupled to the temperature defined at each node, provides a discrete model for simulating 
brittle materials, in agreement with the behavior of continuum elasticity theory. Simulations 
reproduce qualitatively well experiments, but still, a simple description from a theoretical 
point of view is missing. The simulation framework presented here can provide a valuable 
tool to check new proposals for explaining fracture phenomena. Some possible extensions of 
this work follow. 

Here we have analyzed the case in which all the springs have the same force constant, 
or in the elastic language A = \x. This is not the case for most materials. In order to 
simulate materials with different ratios in the Lame coefficients a v^3~ x v^3 reconstruction of 
the triangular lattice with two different force constants may be used |[20|| . 
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It is also tempting to think of this thermally activated cracking effect as a basis of a 
machinery that could microfabricate all kinds of devices. By passing the glass film in one 
direction through the thermal gradient, and then backwards in the other, and also modifying 
the intensity of the temperature change, a large variety of forms can be achieved, and all 
with relatively high accuracy. 

The extension of the calculations done in this work to two dimensional disks is not 
difficult. The question is if this kind of fracture can be accomplished experimentally. 

Another problem of technological interest for which very similar techniques to those 
employed here is stress corrosion cracking. In this case a poisoning agent plays the role of 
temperature by increasing locally the volume of the material (the scalar field that couples 
to elastic displacements is concentration). 
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FIGURES 



FIG. 1. On the left, schematic representation of Yuse and Sano experiment. The right part 
shows the temperature profile in the system (see text). 

FIG. 2. Simplified phase space showing straight, oscillating and several levels of branching for 
systems of size 30x100 and £ = 6. 

FIG. 3. Stresses in the plate, for two different stress conditions that give rise to a straight 
crack: in the upper part thermally driven with AW = 1.8, v = 1, and £ = 6; and in the lower part 
uniform uniaxial stress in the x-axis (<ro = 1.5). Both profiles were taken when the crack is at the 
middle of the plate. System sizes are, for both cases, 30x100. 

FIG. 4. Cross section of the stresses shown in Fig. 3 along the crack. 

FIG. 5. One of the wiggling crack patterns of Fig. 2 (AW = 1.8, and v = 3), and scaling it for 
twice and three times the involved lengths. For twice AW = 1.24, and for three times AW = 1.0. 
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